BNL-NT-08/11 2008/April 



Canonical partition function and finite density phase transition in lattice QCD 

Shinji Ejiri 

Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA 

(Dated: October 10, 2008) 

We discuss the nature of the phase transition for lattice QCD at finite temperature and density. 
We propose a method to calculate the canonical partition function by fixing the total quark number 
introducing approximations allowed in the low density region. An effective potential as a function 
of the quark number density is discussed from the canonical partition function. We analyze data 
obtained in a simulation of two-flavor QCD using p4-improved staggered quarks with bare quark 
mass m/T = 0.4 on a 16^ x 4 lattice. The results suggest that the finite density phase transition at 
low temperature is of first order. 
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I. INTRODUCTION 



The study of the QCD phase diagram at nonzero temperature (T) and quark chemical potential (fiq) is one of 
the most important topics among studies of lattice QCD. In particular, the study of the endpoint of the first order 
phase transition line in the (T, fiq) plane is interesting both from the experimental and theoretical point of view. The 
existence of such a critical point is suggested by phenomenological studies [J, H, Q . The appearance of the critical 
point in the (T, fiq) plane is closely related to hadronic fluctuations in heavy ion collisions and may be experimentally 
examined by an event-by-event analysis of heavy ion collisions. Many trials have been made to find the critical point 
by first principle calculation in lattice QCD 0, H, i, 0, i, i, ES [O, 

Il2l |. However, no definite conclusion on this issue 

is obtained so far. 

One of the interesting approaches is to construct the canonical partition function Zq {T, N) by fixing the total 
quark number (N) or quark number density (p) (isl . 0, [isl . [igI . [TtI . [18| . The canonical partition function is obtained 
from the grand canonical partition function Zgc{T, fj,q) by an inverse Laplace transformation. The relation between 
2Gc{T,Pg) and Zc{T,N) is given by 
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(r,/i,) = / I?[/(dctAf(M,/r))^' e-^« = Yl 2c(r,7V)e^^'/^, (1) 
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where detM is the quark determinant, Sg is the gauge action, and N{ is the number of flavors. iVf in this equation 
must be replaced by N{/4 when one uses a staggered type quark action. In order to investigate the net quark number 
giving the largest contribution to the grand canonical partition function at (T, fj,q), it is worth introducing an effective 
potential VeS as a function of N, 

VMN) ^ ~ lnZc(r, N)^N^^ = - 7V^, (2) 

where / is the Helmholtz free energy. If there is a first order phase transition region, we expect this effective potential 
has minima at more than one value of N. At the minima, the derivative of Vcs satisfies 

Hence, in the first order transition region of T, we expect d{lnZc)/dN{T,N) = —fj.q/T takes the same value at 
different A^. Here, p,*{T,N) is the chemical potential which gives a minimum of the effective potential at {T,N). 

The phase structure in the (T, p) plane and the expected behavior of p* /T are sketched in the left and right panels 
of Fig. [U respectively. The thick lines in the left figure are the phase transition line. We expect that the transition 
is crossover at low density and becomes of first order at high density. Since two states coexist on the first order 
transition line, the phase transition line splits into two lines in the high density region, and the two states are mixed 
in the region between two lines. The expected behavior of p* along the lines A and B are shown in the right figure. 
When the temperature is higher than the temperature at the critical point Tpc (line A), p* increases monotonically 
as the density increases. However, for the case below Tcp (line B), this line crosses the mixed state. Because the two 
states of pi and p2 are realized at the same time, p* does not increase in this region between pi and p2- 

The Glasgow method [3j[lB| has been a well-known method to compute the canonical partition function. Recently, 
such a behavior at a first order phase transition has been observed by Kratochvila and de Forcrand in 4-flavor QCD 
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FIG. 1: Phase structure in the {T,p) plane and the behavior of fiq/T as a function of p. 



with staggered fermions on a 6^ x 4 lattice [l^ calculating the quark determinant by the Glasgow algorithm. However, 
with present day computer resources, the study by the Glasgow method is difficult except on a small lattice. Therefore, 
it is important to consider a method available for a simulation on a large lattice. In this paper, we propose such a 
method for the calculation of the canonical partition function introducing approximations allowed in the low density 
region. 

The method proposed in this paper is based on the following ideas. We adopt a saddle point approximation for 
the inverse Laplace transformation from Zqq to Zq. This approximation is valid when the volume size is sufficiently 
large. We moreover perform a Taylor expansion of hidet M{fiq) in terms of fiq around fiq — and calculate the 
expansion coefficients, as proposed in (T^ . The Taylor expansion coefficients are rather easy to calculate by using the 
random noise method. The saddle point approximation is also based on a Taylor expansion around the saddle point. 
We estimate the Taylor expansion coefficients at the saddle point from the Taylor expansion at /iq = 0. Although we 
must cut off this expansion at an appropriate order in fj,q, this approximation is applicable when the saddle point is 
not far from = 0, and we can estimate the application range where the approximation is valid for each analysis. 

For the calculation of Zc two further technical problems must be solved. The first problem is the problem of 
importance sampling in Monte-Carlo simulations (overlap problem). Since the method proposed in this paper is a 
kind of reweightin g m ethod, the configurations which give important contributions will change when the weight factor 
is changed by fiq ]20l . l2l| . To avoid this problem, we combine configurations generated at many simulation points 
P = Q/g^ covering a wide range of the temperature using a method introduced in (2^ . The second problem is the 
sign problem. We must deal with an expectation value of a complex number in this method. If the fluctuation of the 
complex phase is large, the statistical error becomes larger than the mean value. We use a technique introduced in 
p^ . We consider the probability distribution function in terms of the complex phase of the complex operators. We 
assume the distribution function is well approximated by a Gaussian function and perform the integration over the 
phase. Once this assumption is adopted, the sign problem is completely solved. This assumption is reasonable for 
sufficiently large volume and small fiq/T. 

During the process of this calculation, we will find out why the quark number density changes sharply at the 
transition point and why the density approaches the value of the free quark gas in the high density limit even at low 
temperature. The configurations generated in /i^ = simulations at low temperature are gradually suppressed as the 
density increases. 

In the next section, we explain the method to calculate the canonical partition function using the inverse Laplace 
transformation of Zqq within a saddle point approximation. The problem of the Monte-Carlo sampling is discussed 
in Sec. IIIIl The sign problem is discussed in Sec. IIVI We evaluate d{\n Zc)/dN using data obtained with two-flavors 
of p4-improved staggered quarks in [6]. The result is shown in Sec. |Vl The behavior of d{\n Zc)/dN suggests that 
the phase transition is of first order in the low temperature and high density region. Conclusions are given in Sec. lVIl 
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II. CANONICAL PARTITION FUNCTION 



We calculate the canonical partition function using x Nt lattice and investigate the effective potential V^MN). 
From Eq. ([T]) , the canonical partition function can be obtained by an inverse Laplace transformation [TgI . Ivn . llSj , 
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(4) 



where fj,o is an appropriate real constant and /i/ is a real variable. Note that Zqc{T, fiq + 27riT/3) = Zqc{T, fiq) 
[2^ . The grand canonical partition function can be evaluated by the calculation of the following expectation value at 

= 0. 
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However, with present day computer resources, the exact calculation is difficult except on small lattices. We consider 
an approximation which is valid for large volume and low density. If we select a saddle point as /xg in Eq. ^ when 
the volume is sufficiently large, the information which is needed for the integral is only the value of det M around the 
saddle point. Furthermore, if we restrict ourselves to study the low density region, the value of det M{^q/T) near the 
saddle point can be estimated by a Taylor expansion around //g = 0. The calculations by the Taylor expansion are 
much cheaper than the exact calculations and the studies using large lattices are possible. 

First, we perform the integral in Eq. (jl]) by a saddle point approximation. We denote the quark number density in 
a lattice unit and physical unit as p = N/N^ and p/T^ = pNf , respectively. We assume that a saddle point zq exists 
in the complex Pq/T plane for each configuration, which satisfies 



where (det Af(z)/ det Af (0))^f = exp[7V3D(z)] and D'{z) = dD{z)/dz. 

We then perform a Taylor expansion around the saddle point and obtain the canonical partition function. 



(6) 



Zc{TrpV) = — ZGc(r,o) 



2ti 
3 



2Gc(r,o) 




det AT (0) 



dx 



(T,At,=0) 



exp 



V{D{zo)-pzo--D"{zo)x'' + 



dx 



(T,M,=0) 



--GC 



(r,0)(exp [V (i?(zo)-p2o)]e-'"/2. 
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(7) 



(T,A<,=0) 



Here, D"{z) = d'^D{z)/dz^, V = and D"{z) = \D" {z)\e''°' . We chose a path which passes the saddle point. Higher 
order terms in the expansion of D{z) becomes negligible when the volume V is sufficiently large, since the saddle 
point approximation is a 1/V expansion. 

Next, we calculate the quark determinant by the Taylor expansion around fiq — 0. We define 
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The even derivatives of In det M are real and the odd derivatives are purely imaginary [T9| . The calculation of Z?„ is 
rather easy using the stochastic noise method. D{z), D'{z) and D"{z) in Eq. ^ and O can be evaluated by 



D{z)=NiNtY,Dnz'', D'{z)^NfNtJ2nD„z"-\ D"{z) = NfNtJ2n{n - l)D.nz''^ 



(9) 
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Because Im(Di) ^ Re{D2), the saddle point, i.e. the solution of Eq. is distributed near the real axis and Re(zo) 
increases as p increases. Moreover, for the case that the saddle point zq is on the real axis, the saddle point condition 
is the same as the condition where Re(Z3(z) — pz) is minimized. Hence, exp[yRe(_D(zo) ~ pzo)] in Eq. ([7]) decreases 
exponentially as Re(zo) increases. 
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In this study, we want to focus on the derivative of the effective potential with respect to N or p. Since the effective 
potential VcsiN) is minimized in the thermodynamic limit, i.e. dlogZ^/dN + iig/T = 0, we denote the derivative by 

M; ^ 91nZc(r,jV) ^ l dlnZc{T,pV) 

T dN V dp ' ^ ' 

Within the framework of the saddle point approximation, this quantity can be evaluated by 



zo exp [V (Dizo) - pzo)] e^/^ ^ , 

^ ' U .Mg— (11) 



(exp [V (Dizo) - pzo)] e— /y^TjT^ 



This equation is similar to the formula of the reweighting method for finite pq. The operator in the denominator corre- 
sponds to a reweighting factor, and p^/T is an expectation value of the saddle point calculated with this modification 
factor. We denote the real and imaginary parts of the logarithm of the weight factor by F and 9, 
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F + te = V\ NfNt ^ D^z'S -~pA-\ HV\D"{zo)\] - y . (12) 
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We define the complex phase of the weight factor by 9 and the absolute value of the reweighting factor is exp(i^). 
This weight factor plays an important role around the phase transition point at finite density. 

In the calculation to derive Eq. ([7]), we replaced the order of the path integral of gauge fields and the integral for 
the inverse Laplace transformation. This replacement is essentially important. If one calculates Zq from Zqq using 
Eq. ([4]) after the path integral, the equation which is satisfied at a saddle point is 

N ^pV = d{\nZGc)/d{pjT). (13) 

Hence, in the thermodynamic limit, i.e. when we ignore the finite volume correction, /i* is just equal to the inverse 
function of p{pq) at the saddle point. Therefore, p{pq) must be a discontinuous function or a multivalued function at 
a first order phase transition to obtain the behavior of p*q{p) shown in Fig. [TJ However, if we calculate ln2^GC by a 
Taylor expansion in pq at a temperature in the hadron phase, p{pq) cannot be a discontinuous function. 

In this study, we use Eq. ([TT|) . As we will discuss in detail, we can obtain the behavior of /i*/T suggesting a first 
order phase transition, although the calculations of zq and the weight factor in Eq. (jlip are based on the Taylor 
expansion. The important point is that the weight factor exp(i^ + iO) gives the same effect when the temperature 
changed and configurations which give important contributions to the calculation of Pq/T change gradually as p 
increases. Hence Pq/T does not need to increase monotonously as a function of p even if the saddle point zq is a 
monotonous function of p for each configuration. 



III. MONTE-CARLO ANALYSIS AND REWEIGHTING METHOD FOR /^-DIRECTION 

We analyze data obtained in a simulation of two-flavor QCD using p4-improved staggered quarks with bare quark 
mass m/T = 0.4 @. These data are obtained at 16 simulation points from /3 = 3.52 to 4.00. The corresponding 
temperature normalized by the pseudocritical temperature is in the range ofT/Tc = 0.76 to 1.98, and the pseudocritical 
point {T/Tc = 1) is (ipc ~ 3.65. The ratio of pseudoscalar and vector meson masses is mps/my ~ 0.7 at (i — 3.65. 
The lattice size A^site = -^s x Nt is 16^ x 4. The number of configurations A^conf. is 1000 - 4000 for each (3. Further 
details on the simulation parameters are given in Q. 

We use the data of the Taylor expansion coefficients £)„ up to 0(/i^). The saddle point is found by the following 
procedure. (1) Because absolute values of the odd terms of Z3„ are much smaller than the even terms, we first find 
a solution when the odd terms are neglected, i.e. a solution of NfNt 2nZ)2n^^"~^ — p = 0. The odd terms of £)„ 
are purely imaginary and the even terms are real. Although some fake solutions may appear at large z due to the 
truncation of the higher order terms, the solution is found on the real axis in the low density regime. (2) Next, in the 
vicinity of this solution, we calculate = \NfNt X]„ nDnz"'^ — and find the point where is zero. For the data 
we used in this analysis, the saddle point could be found for every configuration by this procedure except in the very 
low density region of p/T^ < 0.37. Figure [2] shows an example of the distribution of the saddle points, obtained at 
/? = 3.55, 3.63, 3.70 with p/T^ ^ 2.0. 

For the calculation of the derivative of In Zq , the application of the reweighting method for /9-direction is crucial. 
Configurations in a Monte-Carlo simulation are generated with the probability in proportion to the product of the 
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FIG. 2: Distribution of the saddle point at /3 = 3.55, 3.63, 3.70 with p/T^ = 2.0. 
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FIG. 3; Plaquette histogram w{P,0) at /ig = 0. 



weight factor (det M)^'e and the state density of the link fields {U^{x)}. The expectation value is then estimated 
by taking an average of the operator 0[C/;i] over the generated configurations {[/^(x)}. 



(14) 



However, if the value of 0[?7;i] changes very much during a Monte-Carlo simulation and the change of OfC/^] is much 
larger than the size of the probability p{0, /3), the Monte-Carlo method is no longer valid. For example, configurations 
which have large O x p{0,j3) are important for the evaluation of {0)i^p), but such configurations are not generated 
if p{0,/3) is too small. Such a problem occurs in the calculation of Eq. (fTTj) . This problem is called an "overlap 
problem." 

To clarify the problem of the importance sampling, we rewrite Eq. (|lip as 
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P is the plaquette value, 
distribution of P at 



J{exp[F + ie])pw{P,(3)dP 
■)p denotes the expectation value for fixed P dX — 0, and w{P,(3) is the probability 



,(P', /?) = y" VU5{P - P')(det Af(0))^'e-^»('') oc (<5(P - P')) 
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FIG. 4: Expectation value of \zo\ with fixed P for each p/T^. 
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FIG. 5: Expectation vale of F with fixed P for each p/T^. 

This kind of analysis is called the density of state method [l^, HE HI, [13, HI, ■ We define the average plaquette as 
P — —Sg/{6Nsite0) for later discussions. This P is the plaquette value for the standard gauge action, but is a linear 
combination of Wilson loops for improved gauge actions. Because 

_ (XJ(P-PO)(^,^^^o) _ JVUXS{P-P'){detMiO)r^ 
^ (<5(P-PO)(/3.p,=o) !VU5{P-P'){detM{Q)r^ ' ^ ^ 

{X)p is independent of /? for an operator X which does not depend on (3 explicitly. Hence, {X)p can be computed 
at an appropriate p. The probability distribution functions w(P) and —Itiw(P) are given in Fig. [3) We show the 
values of /3 and corresponding T/Tc above these figures. To obtain w{P), we grouped the configurations by the value 
of P into blocks and counted the number of configurations in these blocks. — liiw{P) is normalized by the minimum 
value for each (3. Because the transition from the hadron phase to the quark-gluon phase is a crossover for two-flavor 
QCD, the distribution is always of Gaussian type, and the width of the distribution becomes narrower as the volume 
increases. Moreover, since the suppression factor is exp{6NsitcPP), the peak position of the distribution w(P) moves 
to the right as /3 increases. 

We also calculate the expectation value of \zq\ and F when P is fixed, 

{\zo\)p' - {\zo\SiP P'))/ {6{P P')), {F)p, = {F6{P - P'))/ {5{P ~ P'))- (18) 

The result of (|zo|)p is plotted in Fig. 01 and solid lines in Fig.[5]are (P)p for each p/T^ . (Dashed lines will be explained 
in the next section.) For the calculation of these quantities, we use the delta function approximated by a Gaussian 
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FIG. 6: /3 dependence of — In^oc determined by the consistency condition Eq. (|2ip . 



function, 5{x) w 1/(Ay^) exp[— (a;/A)^]. If making A small, this approximation becomes better but statistical errors 
become larger. Hence, the size of A must be adjusted appropriately. In this study, we adopt A = 0.0025. 

Let us now consider {exp[F + i9])p x w{P) in Eq. Since {F)p increases linearly for P ^ 0.85, Fig. [5] suggests 

that (exp[i^ -I- i9]) p increases exponentially as P increases in this range. Moreover, the slope of (F) p is increasing 
as p/T^ increases. Therefore, for large p/T^, this (exp[i^ + i9])p x w{P) may not decrease even if w{P) decreases 
exponentially at the tail of the distribution generated by a simulation. If the value of (exp[i^ + i9])p x w{P) is still 
large even in the region where the configurations are not generated, the calculation by the Monte-Carlo method is 
completely wrong. Because the width of w{P) becomes narrow for large V, it is essentially important to solve this 
problem if we want to use a large lattice which is required for the saddle point approximation. 

To avoid this problem, we combine all configurations obtained in simulations with many different /3 values, using 
the Ferrenberg and Swendsen method [l^]. The expectation value (O)^ can be calculated from the data obtained by 
more than one simulation point. Pi (i = 1, 2, • • • , Np) by the following equation: 



{OG{p,P)) 



all 



Here, the weight factor G(/3, P) is 



G(/?,P) 



(G(/3,P))an 



(19) 



(20) 



where Nt is the number of configurations at simulation points Pi and (• • • )^y^ means the average over all configurations 
generated at all Pi. The derivation of this equation is given in Appendix A. 

The partition function 2gc(A) is determined by a consistency condition for each i. 



2Gc(A) = (G(ft,P)> 



all 



(21) 



This equation can be solved except for the normalization factor. The result of — ln[ZGc(/3)/2Gc(3.65)] is plotted in 
Fig.H 

We should note that G(/3, P) is independent of the simulation points Pi at which the operators are measured and 
the expectation value is simply given by the average over all configurations generated at many p. If we perform 
simulations at many different P and combine the data, the configurations are distributed in a wide range of P. (See 
Fig. [3l) Among these configurations, important configurations for each calculation are selected by O and G{P,P) 
automatically. This method is particularly important when the volume is large, since the distribution u;(P) is narrow 
if we generate configurations on a large lattice with single p. The overlap problem is solved by this method. However, 
the statistical error is enlarged by the fluctuations of exp[_F + i9] when the density is increased, hence the application 
range of p is determined by the statistical error. Also, we should check that the important configurations arc within 
the range of the plaquette distribution for each calculation. We will discuss this point in Sec. |V] again. 
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IV. SIGN PROBLEM 



Next, we discuss the sign problem that also shows up in the calculation of the canonical partition function. Because 
zq and the additional weight factor in Eq. (fTT|) are complex numbers, the calculation of Eq. (fTTj) suffers from the sign 
problem [U Is^l ■ If the weight factor changes the sign frequently, both the numerator and denominator of Eq. (TTTI) 
become smaller than their statistical errors. To avoid the sign problem, we use a method proposed in [12]. As in 
Eq. (fT2|) , we define the complex phase of the weight factor by 



Im 



V NfNt 



E 

n=l 
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PZQ 



a 
2' 



(22) 



In this definition, is not restricted to the range from —tt to tt because there is no reason that the imaginary part of 
Eq. (|12p must be in the finite range. In fact, this quantity becomes larger as the volume increases. 

It has been discussed in [T^] that histograms of Dn are well approximated by Gaussian functions if a simulation is 
performed at a point away from the critical point with sufficiently large volume. The Taylor expansion coefficients 
in Eq. are given by combinations of traces of products of 9"M/9(/ig/r)" and M~-^. (See the appendix of [Q].) 
Therefore, Z)„ are obtained by the sum of the diagonal elements of such matrices. When the correlation among the 
diagonal elements is small and the volume is sufhciently large, the distribution functions of the expansion coefficients 
and 9 should be of Gaussian type due to the central limit theorem. For example, the diagonal element of the first 
coefficient, lm.[d{\ndet M)/d{iiq/T)] = lm[tT[M~^{dM/d{^q/T))]], is the imaginary part of the local number density 
operator at fiq = 0. If the spatial density correlation is not very strong, a Gaussian distribution is expected. 

Figure[7]is the histogram of the complex phase 9 at p/T^ = 2.0 for the two-flavor QCD simulations with p4-improved 
staggered quarks at /3 = 3.55 (left), 3.63 (middle), and 3.70 (right). These figures suggest that the distribution of 9 is 
well approximated by a Gaussian function. We fit the data of the histogram by a Gaussian function. Dashed lines are 
the fit results. The width of the Gaussian function is different for each distribution obtained by a different parameter. 
If we restrict the phase to the range from — vr to tt by subtracting 27m {n: integer), the complex phase distribution 
is almost flat for the case that the width is much larger than tt, and the flatness indicates the seriousness of the sign 
problem. However, the measurement of the width of the Gaussian distribution is easier than the estimation of the 
flatness of the restricted phase distribution. 

Once we assume a Gaussian distribution for 9, the problem of complex weights can be avoided. We introduce the 
probability distribution as a function of the plaquette P, and 9, 



w{P',F\9') 



VU5{P' - P)5{F' - F)5{9' - 6')(detM(0)) 



(23) 



The distribution function itself is defined as an expectation value at /ig = 0, however F and 9 are functions of p. The 
denominator of Eq. pTjl . (exp[i^ -I- i9]), is given by 



— ( dP I dF ( d9 e^e'%{P,F,9), 

GC J J J 



(24) 



where Zqc — J dP J dF J d9 w{P, F, 9). Because we calculate this expectation value by the reweighting method using 
Eq. p^. the operator in the calculation of Eq. ([M]) is a function of P, F and 9. 
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Since the partition function is real even at nonzero density, the distribution function is symmetric under the change 
from 9 to —9. Therefore, the distribution function is a function of 9^, e.g., il){9) ~ exp[— (02^^ + 04^'' + Og^^ + •••)]. 
And, the distribution function is expected to be well approximated by a Gaussian function when the system size is 
sufficiently large in comparison to the correlation length. We assume the following distribution function in terms of 
6 when P and F are fixed: 



w{P, F,9)^^ ^^iMly^'^p^ F) exp [-02 (P, F)9^] . (25) 
The coefficient 02 {P, F) is given by 

1 . w{P',F\9) 



2a2(P',F') 



J d9 6*2 w{P',F',9) I J 
{9H{P'-P)5{F'-F))^^^^^^^^ _ 



{5{P'-P)5{F'-F))^^^^^^,^ 'iP^n- 
The integration over 9 can be carried out easily and we obtain the denominator of Eq. (jlip 

(^"^")(T,.,^o) - ^/^^/^^/^^ ^'(P,F)e-^%-e- 
— / / w'{P,F)e''e 

■■GC J J 



^{(^')(P.Fr (26) 



/ I?;7e^e-i/(4"2(P.i=^))(detM(0))^'e-^« 
Zgc J 

gFg-l/(4a.(P,F))\ _ (27) 

Since 9 is roughly proportional to the size of the quark matrix M, the value of l/a2 becomes larger as the volume 
increases. Therefore, the phase factor in this quantity decreases exponentially as a function of the volume. However, 
in this framework, the operator in Eq. (|27p is always real and positive for each configuration. This means that the 
expectation value is always larger than its statistical error. Therefore, the sign problem is completely avoided if we 
can assume the Gaussian distribution of 9. The effect from a non-Gaussian term (04) is discussed in Appendix B for 
the case that 04/02 is small. In this case, the effect from 04 is at most 0(/i^), hence the non-Gaussian term may be 
neglected in the low density region even if 04 is nonzero. The complex phase distribution of the quark determinant by 



chiral perturbation theory has been discussed in 31|. The Gaussian distribution is suggested in the 1-loop calculation. 

The complex phase factor for the calculation of (2;oexp[F -|- i9]) in Eq. (fTT|) is calculated repeating the same 
procedure. We consider the complex phase of zo = |zo| exp[i6''], where — tt < 9' < it. Replacing 9 with 9 + 9' in 
Eq. ([77)1 . the suppression factor exp[—{{9 + 9'Y) /2\ is estimated. 

We plot the result of the average of 9^ as a function of P in Fig. [H i.e. 



q2\ 



P' 



H{P-P'))I{5{P-P')). (28) 



We adopt the approximated delta function 5{P) used when we computed (I^oDp and {F)p in Sec. IIIII Although 
we need to average 9'^ as a function of P and F for the calculation of Eq. ([77|) . F dependence is not considered in 
Fig. [H It is found from this figure that {9'^)p decreases linearly as P increases in the range of P < 0.85 and the 
phase fluctuation is small for P > 0.85. This means that the phase factor decreases exponentially as P decreases 
for P < 0.85. This behavior is similar to that of {F)p in Fig. [5] Therefore, the weight factor exp[F — (6'^)(p^i?)/2] 
suppresses the contribution from configurations having small P for large p/T^ ■ 

Moreover, this argument implies that configurations on which the sign problem is serious do not contribute to the 
actual calculations of expectation values. The reason is that the fluctuations of the complex phase {9'^)(^p p-^ are large 
on such configurations and the configurations are suppressed by the weight factor exp[— (^?^)(p j?)/2]. Therefore, even 
if the error due to the Gaussian approximation of the complex phase distribution becomes visible when the phase 
fluctuations are large, the error does not affect to the practical calculations of expectation values so much. 

Here, we should notice that the values of P and F are strongly correlated. We estimate the width of the distribution 
of F for each P and p/T"^ by calculating AF = ^y {{F — {F)p)'^)p. Dashed lines above and below the solid line for 
(F) p in Fig. [5] are the values of (F) p ± AF. Most of the conflgurations characterized by P and F are distributed in 
the narrow region between these two dashed lines. Outside this bound an accurate calculation of {9'^)(^pp^ is difficult, 
since the number of configurations is not enough for the average. However, if we consider that F is approximately 
given as a function of P on each configuration, {9"^) (p^f(p)) will be a function of only P. 
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FIG. 8: Average of the square of the complex phase as a function of P and p/T^ ■ 
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FIG. 9: Expectation value of the plaquette for each pjT^ . Dashed lines are the peak positions of the plaquette distributions 
at pq = 0. 

V. RESULTS AND DISCUSSIONS 

We calculate the slope of InZc, i.e. iJ,*q/T, using Eq. (fTT|) . This quantity is given by the average of the saddle point 
zq multiplying the additional weight factor exp[F + z6']. The volume V — 16"^ is sufficiently large, and we assume that 
the complex phase distribution of the reweighting factor is well approximated by a Gaussian function as discussed in 
Sec IIVI We then replace the weight factor by exp[i^ — (6'^)(i?p)/2] to eliminate the sign problem. We find a saddle 
point zq numerically for each configuration, assuming zq exists near the real axis in the low density region of the 
complex ^q/T plane. 

Before showing the result for i^*/T, it is worth discussing the effect of the weight factor using Eq. ([TS]). The weight 
factor can be approximately estimated by 

(exp[^^ - id])pw{P, 13) « n{P) exp [GN^itePP + {F)p - {0^)p/2] (29) 

because ln(exp[F + i9])p w {F)p — {9'^)p/2 in the leading order and w{P, (3) can be written as fl{P) exp[6j3NsitcP] 
from Eq. (|20p . where ri(P) is the state density in terms of P. 

Since the behavior of {F) p — {9"^) p / 2 for P ^ 0.85 and P ^ 0.85 is different, we consider these two regions separately. 
The configurations below and above P ~ 0.85 are generated in the low and high temperature phases, respectively. In 
the region P > 0.85, the P dependence of {F)p — {0'^)p/2 is small. Therefore, the balance of the weight does not 
change. On the other hand, for P ^ 0.85, {F)p — (6'^)p/2 increases linearly. This has the same effect as when (3 
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changes to 

q2\ 



Peff = P + 77^ 7: J 



P 



\ dP 2 dP J 6K 



(30) 



and the derivatives of (F) p and —{9"^) p increase as p increases. This means that the peak position of the probabihty 
distribution of P, shown in Fig. [3] for fiq = 0, moves to the right as p increases, according to the change of the effective 
(3. This behavior is consistent with our usual expectation, i.e. a phase transition arises when the density is increased 
as well as with increasing temperature. 

For the case that important configurations change, the multiparameter (/?) reweighting in Sec. Illll is effective, since 
the important configurations are automatically selected among all configurations generated at multi-/?, and also this 
method is useful for the interpolation between the simulation points. We combine all data obtained at 16 points of /3 
using the multi-/3 reweighting method. 

In order to observe the change of the important configurations, we calculate the expectation value of the plaquette 
with the additional weight factor in Fig. [51 

{Pexp[F + t0]),^^ 
(exp [i^ + z6'])(T^^^^o) 

The circles in Fig. [9] are (P) at p/T^ — computed without the multi-/3 reweighting method. The solid lines which 
connect these circles are the interpolation using the multi-/3 reweighting method. Of course, (P) for each simulation 
point and that obtained by the multi-/3 reweighting method are consistent with each other at the simulation points. 
However, the line of the interpolation shows waves at some places between the circles. Such a wave appears where 
the configurations are missing in Fig. [31 

For the calculation at finite p/T^, we calculate the phase factor e*^ from {6'^)p in Fig. [51 The F dependence of 
{9'^){p,F) is neglected because the values of F is approximately given as a function of P as seen in Fig. [51 The results 
for (P) at p/T^ = 2.0,4.0 and 6.0 are plotted from below. This quantity indicates the plaquette value of the most 
important configuration when the weight factor is modified. 

From this figure, we find that (P) becomes larger for (3 < fipc = 3.65 and does not change very much for /3 > 3.65 
when p/T^ is increased. This means that, even when we measure an expectation value at small temperature (small (3), 
the configurations generated by simulations at higher temperature (larger P) are used for the measurement at finite 
p/T^ ■ Moreover, for the measurement at sufficiently large p/T^ , the configurations generated in the low temperature 
phase are completely suppressed by the additional weight factor even when the temperature is small. This property 
explains why the phase transition happens when the density is increased while keeping the temperature fixed. 

However, the errors on these results become larger as p/T^ increases. One of the reasons may be the missing 
configurations between the peaks of the plaquette distributions. As seen in Fig. [31 the configurations are not distributed 
uniformly in the range of P which is necessary in this analysis, and correct results cannot be obtained if the important 
configurations are missing. At low temperature, the important value of P changes very much as p increases, therefore 
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we calculate ^*/T only when the expectation value of P is at the peak positions of the plaquette distributions in 
Fig. [31 Dashed lines in Fig [5] are the peak positions. 

Another important point is that (|zo|)p and {F) p are strongly correlated with each other. Figure [4] is the average 
of |zo| as a function of the plaquette value of each configuration. \zq\ increases as P decreases. Because Di is purely 
imaginary, —F = Re[V{pzQ — N[NtDiZo)] + 0{zq) becomes large as Re(zo) increases, which is seen in Fig. [5l and 
the contribution from the configurations which have large zq is suppressed by the additional weight factor. Although 
the value of |zo| for each configuration increases monotonically as a function of p/T^, nontrivial behavior in fi*/T is 
expected due to the suppression factor. 

We plot the result of /z*/T in Fig. [10] as a function of p/T^ for each temperature T/Tc{l3). Dashed lines are cubic 
spline interpolations of these results. The dot-dashed line is the value of the free quark-gluon gas in the continuum 
theory. 



(32) 



In this calculation, we neglected the F dependence of {0'^)(p^f) because the values of P and F are strongly correlated. 
The systematic error due to this approximation is discussed in Appendix C. The error seems to be small. 

From Fig. [TOl we find that a qualitative feature of changes around T/Tc ~ 0.8, i.e. /i*/T increases mono- 

tonically as p increases above 0.8, whereas it shows an s-shape below 0.8. This means that there is more than one 
value of p/T^ for one value of /i^/T below T/Tc ~ 0.8. This is a signature of a first order phase transition. The 
critical point in the (T, pq) plane is estimated in [l2| by calculating the effective potential in terms of the plaquette 
value using the same configurations. The estimation is {T/Tc, Pq/T) ~ (0.76,2.5). Because the estimation from the 
effective potential is rather ambiguous, the difference between the new and old results of the critical temperature may 
be a systematic error. The critical value of Pq/T is about 2.4. This is almost consistent with the previous result by 
the different method. The error from the truncation of the Taylor expansion of the quark determinant is discussed 
in [l2|. The difference between the results at 0{pq) and 0{pq) is found to be small at Pq/T — 2.5 for the data 
used in this study. Therefore, the error due to the truncation would not affect the qualitative conclusions. Although 
further studies including justifications of the approximations used in this analysis are necessary for more quantitative 
investigations, this result suggests the existence of the first order phase transition line in the {T,pq) plane. 



VI. CONCLUSIONS 



We studied the canonical partition function as a function of p/T^ performing an inverse Laplace transformation. 
We analyzed the data obtained with two-flavors of p4-improved staggered quarks in Q and calculated the derivative 
of the canonical partition function with respect to p. The problems in this calculation were discussed. To avoid the 
problems, we adopted the following approximations. First, we estimate the quark determinant from the data of a 
Taylor expansion up to 0{pq) because the direct calculation of the quark determinant is still difficult except on a 
small lattice. Although terms of higher than p^ are omitted, this analysis is valid in the low density region. Second, 
we use a saddle point approximation for the inverse transformation, assuming the volume is sufficiently large. Third, 
we assume that the probability distribution of the complex phase of the operator in the calculation of p* /T^ can be 
well approximated by a Gaussian function. 

Using multiparameter reweighting method, we combined the configurations generated by pq — simulations at 16 
simulation points (/3) which cover a wide range of the temperature. It is found that the increase of Pq/T^ becomes 
larger as the temperature decreases in the low density region. However, the contribution from the configurations 
generated at low temperature gradually decreases in the measurement of Pq/T^ as p/T'^ increases even at low temper- 
ature. And, Pq/T^ approaches the value of the free quark gas in the high density limit for all temperatures investigated 
in this study. The most interesting result is that Pq/T^ as a function of p/T^ shows an s-shape at T ^ 0.8. This 
means that the effective potential Kff in terms of the density has two minima. Therefore, this result strongly suggests 
the existence of the first order phase transition line in the low temperature and high density region. Since the data 
we used in this study is obtained by a simulation with much heavier quark masses than the physical quark masses, 
simulations near the physical mass point are very important. It is also necessary to increase the accuracy of the 
approximations we have used in this study. 
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Appendix A: Multi-parameter (/3) reweighting method 

We discuss a method to combine configurations obtained by simulations with many /3, following the method by 
Ferrenberg and Swendsen p^ . We define a /3-independent distribution function of the plaquette value, 

n{P')^ I VU6{P - P'){detM{Q)f'. (33) 



The relation between Vl{P) and w{P, f3) in Eq. (fT6| is n{P) cxp{6f3NsitcP) = w{P,(3). Then, the expectation value of 
an operator C as a function of P is given by the equation 

(O), - ^ / 0(P)17(P)e6/'A'-.Prfp ^ 1 0^ = I n(P)e6/3A^=-^dP, (34) 

where A^conf . is the number of configurations, and {conf } denotes the sum of O over all configurations generated 
in a simulation at /3. 

The expectation value (O) ^ is also calculated from the data obtained by more than one simulation point, (3i 
(i = 1, 2, • • • , N/3). The Eq. ^ is evaluated by 

/ o(p)r!(p)e6'3^='-^dP = / 0(p)r>(p) dP 

i=l /3i,{conf.} 2^j=\^^3^ ^GCyPj) 

where Ni is the number of configurations at simulation points (3i. Hence, 



F y (35) 



^ (G(/3,P)).n ' ^ ^ 



Here, the weight factor G(/3, P) is 



G(/3,P) = ^v ' (37) 

and (• • • )aU means the average over all configurations generated at all Pi. 

The partition function ZQciPi) is determined by a consistency condition for each i, 

2Gc(ft) = E E G(A,P) = (G(A,P))aii. (38) 

J = l /3j,{conf.} 

This equation can be solved except for the normalization factor. We should note that G{(3, P) is independent of 
the simulation points (3i at which the operators are measured, and important configurations for each calculation are 
selected by the weight factor automatically. 

Appendix B: Effect from non-Gaussian terms in the phase factor 

We estimate the phase factor when the distribution is slightly different from Gaussian. We consider a distribution 
function with small a4(P, P)/a2(P, P), 



-1 

0-2 I 3a4 



, 2 

04 



0.2 



i5'(P,P)e-("=^'+'^*^*^ (39) 
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In this case, the phase factor, exp[— l/(4a2)], changes to 



1 

4a2 



(40) 



and also the expectation value of 9^ for fixed P and F becomes 



1 



3a4 



){P.F) 



2a2 2a% 



O 



(41) 



From this equation, the term of 804/(402) in Eq. (PD)) is absorbed into (0^)/2, hence the leading contribution from 04 
in the phase factor is exp[— 04/(16o|)]. The value of 04 can be evaluated by the Binder cumulant, 



iP.F) 



^^)(P,F) 



(42) 



Because On ^ 



0{fj,g) for the chemical potential at the saddle point zo, the effect from 04 becomes larger as 



the density increases. Therefore, for the case of 04/02 ^ 0(1) in the ^ limit, exp[— 04/(160!)] is 0{p^) at most. 



I.e. 



= (e^exp[-(02)(P,;.)/2 + O(/i6)]) 



(T,p<,=0) 



(43) 



This argument suggests that the approximation by the Gaussian distribution is valid for the investigation of the low 
density region even if 04 is nonzero, however the estimation of the range of p in which the non-Gaussian contribution 
is small may be important as well as the application range of the Taylor expansion in /ig. 



Appendix C: Suppression factor from complex phase fluctuation 

In the calculation of fJ-^/T, we need to calculate the suppression factor from the complex phase fluctuation, 
exp[—{9'^){P,F,p/T^)]. This factor should be a function of P, F and p/T^, however we neglected the F depen- 
dence in the calculation of Fig. 1101 In this section, we discuss the error from this approximation. As shown in Fig. O 
the values of P and F are strongly correlated. The solid line is the mean value of F among the configurations having 
the plaquette value P for each p/T^. The dashed lines above and below the solid line show the fluctuations. Most 
of the configurations characterized by P and F are distributed in the narrow region between the two dashed lines. 
Therefore, we have approximated (0^) as a function of only P and p/T"^ in Fig. [51 

To estimate the importance of the F dependence, we calculate Pq/T using another approximation of {9'^) which 
changes according to F of each configuration. From the data of (F) p as a function of p/T^ for each P, which is shown 
in Fig. [51 we find p/T^ which gives F for each configuration and find the value of {9'^) at this p/T^ using the data 
of Fig. [51 We then obtain (6*^) as a function of P and F, but the p/T^ dependence is neglected. We calculate Pg/T 
using this {9'^){p,f) and compare the previous result to estimate the systematic error due to the approximation in 
(0^). The result is shown in Fig. [Til The dotted lines are the spline interpolation in Fig.[Tni The difference between 
the results of Pg/T in Figs. [Till and [TT] seems to be small. This suggests that the determination of (6*^) with three 
parameter {P^ F, p/T^) is not very important for the qualitative argument. 

Note: In the jackknife error estimation of this calculation, we have neglected the dispersion of {9'^)(p^f) among the 
jackknife ensemble. Therefore, statistical errors in Fig. [11] are smaller than those in Fig. [TOl The small error does not 
mean that the analysis in this appendix gives better results with smaller statistical error. In fact, the errors in Fig. 1101 
become the same size if we neglect the dispersion of (6'^) in the jackknife analysis. 
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